import numpy as np
import matplotlib.pyplot as plt

class Atmosphere(object):
    g = 9.8
    P0 = 101300.0
    T0 = 288.16
    rho0 = 1.225
    Mu0 = 1.825e-006
    Gamma = 1.4
    Rs = 8314.32
    M0 = 28.9644
    Density = [1.0, -9.59387E-005, 3.53118E-009, -5.83556E-014, 2.28719E-019]
    Pressure = [1.0, -0.000118441, 5.6763E-009, -1.3738E-013, 1.60373E-018]
    Temperature = [1.0,-2.27712e-005, 2.18069E-010, -5.71104E-014, 3.97306E-018]
        
    def __init__(self):
        pass
    
    @staticmethod
    def polyn(af, f):
        return (((af[4]*f+af[3])*f+af[2])*f+af[1])*f+af[0]
    
    @staticmethod
    def pressure(h):
        if(h > 18300.0):
            return (18300.0/h)*P0*Atmosphere.polyn(Atmosphere.Pressure,18300.0)
        else:
            return Atmosphere.P0*Atmosphere.polyn(Atmosphere.Pressure,h)
    
    @staticmethod
    def temperature(h):
        if(h > 18300.0):
            h = 18300.0
        return Atmosphere.T0 * Atmosphere.polyn(Atmosphere.Temperature, h)
    
    @staticmethod
    def density(h):
        if(h > 18300.0):
            return (18300.0 / h) * Atmosphere.rho0 * Atmosphere.polyn(Atmosphere.Density, 18300.0)
        else:
            return Atmosphere.rho0 * Atmosphere.polyn(Atmosphere.Density, h)
    
    @staticmethod    
    def viscosity(h):
        return Atmosphere.Mu0 * pow(Atmosphere.temperature(h)/Atmosphere.T0, 0.76000000000001)
    
    @staticmethod
    def kineticViscosity(h):
        return Atmosphere.viscosity(h) / Atmosphere.density(h)
    
    @staticmethod
    def speedOfSound(h):
        t = Atmosphere.temperature(h)
        return np.sqrt((Atmosphere.Gamma*Atmosphere.Rs/Atmosphere.M0)*t)


if __name__ == '__main__':
    print Atmosphere.g
    print Atmosphere.P0
    print Atmosphere.T0
    print Atmosphere.rho0
    print Atmosphere.Mu0
    print Atmosphere.pressure(0)
    print Atmosphere.temperature(0)
    print Atmosphere.density(0)
    print Atmosphere.viscosity(0)
    print Atmosphere.kineticViscosity(0)
    print Atmosphere.speedOfSound(0)
    
    alts = np.linspace(0,18300,1000)
    sigma = []
    omega = []
    for alt in alts:
        sigma.append(Atmosphere.density(alt)/Atmosphere.rho0)
        omega.append(Atmosphere.pressure(alt)/Atmosphere.P0)
    plt.plot(sigma,alts)
    plt.plot(omega,alts)
    plt.grid()
    plt.show()
    
